/**
 * Implement the vector computation. 
 * The classes Vector and Matrix are provided. 
 * @author Songgang Xu (sxu at tamu dot edu)
 * @file matrix.h
 */

#ifndef XMESH_VECTOR_H
#define XMESH_VECTOR_H

#include <algorithm>
#include <iostream>
#include <iterator>
#include <functional>
#include <cmath>
#include <numeric>
#include <vector>
#include <memory>

using namespace std ;


/**
 * Define Vector.  
 */
template <typename ValueT>
class Vector
{
protected:
	typedef ValueT ValueType ;
	typedef Vector<ValueT> VectorType ;
public:
	typedef Vector<ValueT> PointType ;
public:
	/**
	 * Constructor. 
	 */
	Vector ():m_size(0) {}

	/**
	 * Copy Constructor. 
	 * @param _v Vector instance. 
	 */
	Vector (const Vector & _v): m_size(_v.m_vector.size()), m_vector (_v.m_vector)  {}

	/**
	 * Constructor. 
	 * @param _s Size of vector. 
	 */
	Vector (unsigned _s): m_size(_s) {m_vector.resize(m_size); }

	/**
	 * Constructor for 2D vector. 
	 *
	 */
	Vector (const ValueType & _v1, const ValueType & _v2) : m_size(2)
	{
		m_vector.resize(m_size);
		m_vector[0] = _v1 ;  m_vector[1] = _v2 ; 
	}

    /**
	 * Constructor for 3D vector. 
	 *
	 */
	Vector (const ValueType & _v1, const ValueType & _v2, const ValueType & _v3) : m_size(3)
	{
		m_vector.resize(m_size);
		m_vector[0] = _v1 ;  m_vector[1] = _v2 ; m_vector[2] = _v3 ;
	}

	/**
	 * Constructor for nD vector. 
	 */
	Vector ( unsigned _n , const ValueType * _v) : m_size (_n)
	{
		m_vector.resize(m_size) ;
		for (unsigned i = 0; i < m_size; ++i)
		{
			m_vector[i] = _v[i]; 
		}
	}

    /**
	 * Get vector size. 
	 */
	unsigned size() const { return this->m_size;}

	/**
	 * Resize the vector.
	 * @param _s
	 */
	void resize (unsigned _s)
	{
		this->m_vector.resize(_s); 
		this->m_size = _s ;
	} 
	 
	/**
	 * Clear the vertor. 
	 */
	void clear ()
	{
		this->m_vector.clear();
		this->m_size = 0; 
	}

	/**
	 * Reference of vector. 
	 * @param _i Index of element. 
	 */
	ValueT & operator [] (unsigned _i) {return this->m_vector[_i];}

	/**
	 * Value of element in vector. 
	 * @param _i Index of element. 
	 */
	const ValueT & operator [] (unsigned _i) const {return this->m_vector[_i];}

	/**
	 * Overload operator = (ValueT) for vector initialization. 
	 * @param _v Value. 
	 * @return Vector reference. 
	 */
	Vector & operator = (const ValueT & _v) 
	{
//		this->m_idxcount = 0; 
		this->m_vector.clear(); 
		this->m_vector.push_back( _v );
		this->m_size = 1 ;
		return (*this); 
	}

	/**
	 * Overload operator comma  for vector initialization.
	 * @param _v Value.
	 * @return Vector reference. 
	 */
	Vector & operator , (const ValueT & _v)
	{
//		++this->m_idxcount ;
//		if (this->m_idxcount >= this->m_size) throw ("Vector initialization exceeds bound. ");
//		this->m_vector[m_idxcount] = _v; 
		this->m_vector.push_back( _v );
		++ (this->m_size ) ;
		return (*this);
	}

	/**
	 * Append a value into the vector.
	 * @param _v Value. 
	 */
	void append (const ValueT & _v)
	{
		this->m_vector.push_back(_v);
		++ (this->m_size) ;
	}

	/**
	 * Make every component to be non-negative. 
	 * @return None.
	 */
	void positive ()
	{
		for (unsigned i = 0; i < m_size; ++i)
		{
			m_vector[i] = abs (m_vector[i]);
		}
	}
public:

	/**
	 * Calculate the module square of Vector. 
	 */
	ValueT square_module ()
	{
		ValueT sum = ValueT(0);
		for (unsigned i = 0; i < m_size ; ++i)
		{
			sum += this->m_vector[i] * this->m_vector[i]; 
		}
		return sum ;
	}

    /**
	 * Calculate the module of vector. 
	 */
	ValueT module ()
	{
		return sqrt (square_module());
	}

	/**
	 * Normalize the vector. 
	 */
	Vector & normalize ()
	{
		ValueT m = module (); 
		for (unsigned i = 0; i < m_size; ++i)
		{
			this->m_vector[i] /= m ;
		}
		return (*this);
	}
	


    /**
	 * Overload of operator *
	 * @param lhs Left hand operand. 
	 * @param rhs Right hand operand. 
	 */
	friend inline ValueT operator * (const Vector & lhs, const Vector & rhs)  
	{
		if (!(lhs.size() == rhs.size())) throw ("The size of vectors doesn't match. ");

		ValueT result = ValueT(0);
		for (unsigned i = 0 ;i < lhs.size(); ++i)
		{
			result += lhs[i] * rhs[i];
		}
		return result; 
	}

    /**
	 * Overload of operator +=
	 * @param _v Vector. 
	 * @return incremented Vector. 
	 */
	inline Vector & operator += (const Vector & _v)
	{
		if (!(this->size() == _v.size())) throw ("The size of vectors doesn't match. ");
		for (unsigned i = 0 ;i < this->size(); ++i)
		{
			 this->m_vector[i] += _v[i];
		}
		return (*this); 
	}

    /**
	 * Overload of operator +
	 * @param lhs Left hand operand. 
	 * @param rhs Right hand operand. 
	 */
	friend inline Vector operator + (const Vector & lhs, const Vector & rhs) 
	{
		if (!(lhs.size() == rhs.size())) throw ("The size of vectors doesn't match. ");
		Vector result (lhs);
		result += rhs; 
		return result; 
	}

    /**
	 * Overload of operator -
	 * @param lhs Left hand operand. 
	 * @param rhs Right hand operand. 
	 */
	friend inline Vector operator - (const Vector & lhs, const Vector & rhs) 
	{
		if (!(lhs.size() == rhs.size())) throw ("The size of vectors doesn't match. ");
		Vector result (lhs.size());
		for (unsigned i = 0 ;i < lhs.size(); ++i)
		{
			result [i] = lhs[i] - rhs[i];
		}
		return result; 
	}

    /**
	 * Overload of operator *
	 * @param lhs Left hand operand. 
	 * @param _v Scale value. 
	 */
	friend inline Vector operator * (const Vector & lhs, const ValueT & _v ) 
	{
		Vector result (lhs.size());
		for (unsigned i = 0 ;i < lhs.size(); ++i)
		{
			result [i] = lhs[i] * _v;
		}
		return result; 
	}

    /**
	 * Overload of operator *
	 * @param _v Scale value. 
	 * @param rhs Left hand operand. 
	 */
	friend inline Vector operator * ( const ValueT & _v, const Vector & rhs )
	{
		Vector result (rhs.size());
		for (unsigned i = 0 ;i < rhs.size(); ++i)
		{
			result [i] = rhs[i] * _v;
		}
		return result; 
	}

    /**
	 * Overload of operator /=
	 * @param _v Vector. 
	 * @return incremented Vector. 
	 */
	inline Vector & operator /= (const ValueT & _v)
	{
		for (unsigned i = 0 ;i < this->size(); ++i)
		{
			 this->m_vector[i] /= _v;
		}
		return (*this); 
	}

    /**
	 * Overload of operator /
	 * @param lhs Left hand operand. 
	 * @param _v Scale value. 
	 */
	friend inline Vector operator / (const Vector & lhs, const ValueT & _v ) 
	{
		Vector result (lhs);
		result /= _v; 
		return result; 
	}

	/**
	 * Overload of operator ^. Cross product of two vectors. 
	 * @param lhs Left hand operand. 3D vector.
	 * @param rhs right hand operand. 3D vector. 
	 */
	friend inline Vector operator ^ (const Vector & lhs, const Vector & rhs) 
	{
		if (!(lhs.size() == 3 && rhs.size() == 3)) throw ("Cross product is computed from two 3D vectors. ");
		return Vector (lhs[1] * rhs[2] - lhs[2] * rhs[1], 
					   lhs[2] * rhs[0] - lhs[0] * rhs[2],
					   lhs[0] * rhs[1] - lhs[1] * rhs[0] );
	}

    /**
	 * Overload of ostream operator
	 * @param out Output stream.  
	 * @param _v Vector. 
	 */
	friend inline std::ostream& operator<<(std::ostream & out, const Vector & _v) 
	{
		out <<'[';
		for (unsigned i = 0; i < _v.m_size - 1; ++i) out << _v.m_vector[i] <<',';
		if (_v.m_size > 0) out <<_v.m_vector[_v.m_size -1]<<']';
		return out;
	}

    /**
	 * Overload of istream operator
	 * @param in Input stream.  
	 * @param _v Vector. 
	 */
	friend inline std::istream& operator>>(std::istream& in, Vector& _v) 
	{
		for (unsigned i = 0; i < _v.size(); ++i)
		{
			in >> _v.m_vector[i]; 
		}
		return in ;
	}

	/**
	 * Overload comparison operator <.
	 * We order two vectors having the same dimension in the alphabetic order. 
	 */
	friend inline bool operator < (const Vector & lhs, const Vector & rhs)
	{
		for (unsigned i = 0; i < lhs.size(); ++i)
		{
			if (lhs[i] < rhs[i]) return true; 
			if (lhs[i] > rhs[i]) return false; 
		}
		return false; 
	}
	
	/**
	 * Overload comparison operator >.
	 * We order two vectors having the same dimension in the alphabetic order. 
	 */
	friend inline bool operator > (const Vector & lhs, const Vector & rhs)
	{
		return  (rhs < lhs) ;
	}

	/**
	 * Overload of operator == 
	 * @param lhs Left hand operand. 
	 * @param rhs Right hand operand. 
	 */
	friend inline bool operator == (const Vector & lhs, const Vector & rhs)  
	{
		return (lhs.m_vector == rhs.m_vector);
	}

	/**
	 * Overload of operator != 
	 * @param lhs Left hand operand. 
	 * @param rhs Right hand operand. 
	 */
	friend inline bool operator != (const Vector & lhs, const Vector & rhs)  
	{
		return (!(lhs == rhs));
	}

	/**
	 * Overload comparison operator <=.
	 * We order two vectors having the same dimension in the alphabetic order. 
	 */
	friend inline bool operator <= (const Vector & lhs, const Vector & rhs)
	{
		return  !(rhs < lhs) ;
	}

	/**
	 * Overload comparison operator >=.
	 * We order two vectors having the same dimension in the alphabetic order. 
	 */
	friend inline bool operator >= (const Vector & lhs, const Vector & rhs)
	{
		return  !(lhs < rhs) ;
	}
private:
	unsigned m_size; 
	std::vector<ValueT> m_vector; 
	
	unsigned m_idxcount; 
};


/**
 * Define a 2D matrix class. 
 */
template <typename ValueT>
class DenseMatrix
{
protected:
	/** 
	 * Type name for value 
	 */
	typedef ValueT value_type ;
	typedef Vector<ValueT> VectorType; 

public:

	/**
	 * Constructs a matrix class. 
	 */
	DenseMatrix (): m_row (0), m_col(0){}

	/**
	 * Constructs a matrix class of r x c.
	 * @param _r row of the matrix. 
	 * @param _c col of the matrix. 
	 */
	DenseMatrix (unsigned _r, unsigned _c) : m_row(_r), m_col (_c) 
	{
		allocate_memory();
	}

	/**
	 * Constructs a matrix class of r x c with the initial value v.
	 * @param _r row of the matrix. 
	 * @param _c col of the matrix. 
	 * @param _v initial value. 
	 */
	DenseMatrix (int _r, int _c,  const ValueT  & _v) : m_row(_r), m_col (_c) 
	{
		allocate_memory (this->m_pElements, m_row, m_col, _v);
	}

	/**
	 * Constructs a matrix from a vector. The vector is viewed to be column vector. 
	 */
	DenseMatrix (const VectorType & _v): m_row (_v.size()), m_col (1)
	{
		allocate_memory();
		for (unsigned i = 0; i < m_row; ++i)
		{
			m_pElements[i][0] = _v[i];
		}
	}

	/**
	 * Copy Constructor.
	 * @param _m DenseMatrix Object. 
	 */
	DenseMatrix (const DenseMatrix & _m)
	{
		m_row = _m.get_row();
		m_col = _m.get_col();
		allocate_memory(); 
		for (unsigned i = 0; i < m_row; ++i)
		{
			for (unsigned j = 0; j < m_col; ++j)
			{
				m_pElements[i][j] = _m.m_pElements[i][j];
			}
		}
	}

	/**
	 * Overload of operator = 
	 * @param _m DenseMatrix Object. 
	 */

	DenseMatrix & operator = (const DenseMatrix & _m)
	{
		//if (this->m_pElements != NULL) deallocate_memory();
		m_row = _m.get_row();
		m_col = _m.get_col();
		allocate_memory(); 
		for (unsigned i = 0; i < m_row; ++i)
		{
			for (unsigned j = 0; j < m_col; ++j)
			{
				m_pElements[i][j] = _m.m_pElements[i][j];
			}
		}
		return (*this) ;
	}

	/**
	 * Destructor. 
	 */
	virtual ~DenseMatrix () {deallocate_memory(); }


	/**
	 * Resize matrix. 
	 * @param _row New row value.
	 * @param _col New col value. 
	 * @return None. 
	 */
	void resize (unsigned _row, unsigned _col)
	{
		/**Fix me here*/
		this->m_row = _row;
		this->m_col = _col; 
		this->m_pElements.resize(_row);
		for (unsigned i = 0; i < _row; ++i)
		{
			m_pElements[i].resize(_col);
		}
	}

    /**
	 * Overload operator [].
	 * @param _i Index of the first dimension
	 * @return Reference of row vector. 
	 */
	Vector<ValueT> & operator [] (unsigned _i)
	{
		return m_pElements[_i];
	}

    /**
	 * Overload operator [].
	 * @param _i Index of the first dimension
	 * @return Constant reference of row vector. 
	 */
	const Vector<ValueT> & operator [] (unsigned _i) const 
	{
		return m_pElements[_i];
	} 


	/**
	 * Overload operator () to fetch element as L-Value. 
	 * @param _r Row value.
	 * @param _c Col value.
	 * @return Reference of element. 
	 */
	ValueT & operator () (unsigned _r, unsigned _c)
	{
		return m_pElements[_r][_c];
	} 
    /**
	 * Overload operator () to fetch element. 
	 * @param _r Row value.
	 * @param _c Col value.
	 * @return Constant reference of element. 
	 */
	const ValueT & operator () (unsigned _r, unsigned _c) const
	{
		return m_pElements[_r][_c];
	}
	/**
	 * Get the row of matrix. 
	 * @return the row count of matrix. 
	 */
	unsigned get_row () const {return this->m_row; }

	/**
	 * Get the col of matrix. 
	 * @return the col count of matrix. 
	 */
	unsigned get_col () const {return this->m_col; }

	/**
	 * Set the row of matrix. 
	 * @param _r Row value. 
	 */
	void set_row (unsigned _r ) {this->m_row = _r; }

	/**
	 * Set the col of matrix.
	 * @param _c Col value. 
	 */
	void set_col (unsigned _c ) {this->m_col = _c; }
	
	//unsigned get_dim1 () const {return this->m_row; }
	//unsigned get_dim2 () const {return this->m_col; }

	/**
	 * Check if the matrix is empty. 
	 * @return if the matrix is empty. 
	 */
	bool is_empty() const 
	{
		return (m_pElements == NULL && m_row == 0 && m_col == 0) ;
	}

	/**
	 * Check if the matrix is square.
	 * @return if the matrix is square. 
	 */
	bool is_square () const 
	{
		return (m_row == m_col);
	}

	/**
	 * Check if two matrices have the same dimension. 
	 * @return if two matrices have the same dimension.
	 */
	friend bool is_dim_match(const DenseMatrix & _m1, const DenseMatrix & _m2)
	{
		return (_m1.get_row () == _m2.get_row() && _m1.get_col () == _m2.get_col());
	}

public:
	/**
	 * Overload plus operator. 
	 * @param lhs Left hand operand.
	 * @param rhs Right hand operand. 
	 * @return the sum of two matrix. 
	 */
	friend DenseMatrix operator + (const DenseMatrix & lhs, const DenseMatrix & rhs)
	{
		if (! (lhs.get_row() == rhs.get_row() && lhs.get_col() == rhs.get_col() ))
		{
			throw ("Matrices dimensions do not match.");
		}
		unsigned _row = lhs.get_row(); 
		unsigned _col = lhs.get_col ();
		DenseMatrix result(_row, _col); 
		for (int i = 0; i < _row; ++i)
		{
			for (int j = 0; j < _col; ++j)
			{
				result (i,j ) = lhs(i,j ) + rhs(i,j);
			}
		}
		return result ;
	}

	/**
	 * Overload minus operator. 
	 * @param lhs Left hand operand.
	 * @param rhs Right hand operand. 
	 * @return the difference of two matrix. 
	 */
	friend DenseMatrix operator - (const DenseMatrix & lhs, const DenseMatrix & rhs) 
	{
		if (! (lhs.get_row() == rhs.get_row() && lhs.get_col() == rhs.get_col() ))
		{
			throw ("Matrices dimensions do not match.");
		}
		unsigned _row = lhs.get_row(); 
		unsigned _col = lhs.get_col ();
		DenseMatrix result(_row, _col); 
		for (int i = 0; i < _row; ++i)
		{
			for (int j = 0; j < _col; ++j)
			{
				result (i,j ) = lhs(i,j ) - rhs(i,j);
			}
		}
		return result ;
	}

	/**
	 * Overload times operator. 
	 * @param lhs Left hand operand.
	 * @param rhs Scale value. 
	 * @return the product of the matrix and the scale value. 
	 */
	friend DenseMatrix operator * (const DenseMatrix & lhs, ValueT & _v) 
	{
		unsigned _row = lhs.get_row(); 
		unsigned _col = lhs.get_col ();
		DenseMatrix result(_row, _col); 
		for (int i = 0; i < _row; ++i)
		{
			for (int j = 0; j < _col; ++j)
			{
				result (i,j ) = lhs(i,j ) * _v;
			}
		}

		return result ;
	}

	/**
	 * Overload times operator. 
	 * @param lhs Scale value. 
	 * @param rhs Right hand operand. 
	 * @return the product of the matrix and the scale value. 
	 */
	friend DenseMatrix operator * (ValueT & _v, const DenseMatrix & rhs) 
	{
		unsigned _row = rhs.get_row(); 
		unsigned _col = rhs.get_col ();
		DenseMatrix result(_row, _col); 
		for (int i = 0; i < _row; ++i)
		{
			for (int j = 0; j < _col; ++j)
			{
				result (i,j ) = rhs(i,j ) * _v;
			}
		}
		return result ;
	}


	/**
	 * Overload times operator. 
	 * @param lhs Matrix.
	 * @param rhs Vector. 
	 * @return the product of the matrix and the vector. 
	 */
	friend DenseMatrix operator * (const DenseMatrix & lhs, const VectorType & _v)
	{
		unsigned _row = lhs.get_row(); 
		unsigned _col = lhs.get_col ();
		if ( !(lhs.get_col() == _v.size())) 
		{
			throw ("Dimensions of matrix and vector do not match. ");
		}
		DenseMatrix vector_matrix (_v);

		return (lhs * vector_matrix);
	}

	/**
	 * Overload times operator. 
	 * @param lhs Left hand operand.
	 * @param rhs Right hand operand. 
	 * @return the product of two matrix. 
	 */
	friend DenseMatrix operator * (const DenseMatrix & lhs, const DenseMatrix & rhs) 
	{
		if (! (lhs.get_col() == rhs.get_row() ))
		{
			throw ("Matrices dimensions do not match.");
		}

		unsigned _row = lhs.get_row(); 
		unsigned _col = rhs.get_col ();
		unsigned _k = lhs.get_col ();
		DenseMatrix result(_row, _col); 
		for (unsigned i = 0; i < _row; ++i)
		{
			for (unsigned j = 0; j < _col; ++j)
			{
				result(i,j) = 0; 
				for (unsigned k = 0; k < _k ; ++k)
				{
					result(i,j) += lhs(i,k ) * rhs(k,j);
				}
			}
		}
		return result ;
	}
	
	/**
	 * Transpose the matrix.
	 * @param Matrix to be transposed.
	 * @return Transposed matrix. 
	 */
	friend DenseMatrix transpose (const DenseMatrix & _m) 
	{
		unsigned _col = _m.get_row(); 
		unsigned _row = _m.get_col ();
		DenseMatrix result(_row, _col); 
		for (int i = 0; i < _row; ++i)
		{
			for (int j = 0; j < _col; ++j)
			{
				result (j, i) = _m(i, j) ;
			}
		}
		return result ;
	}


public: 

	/**
	 * Output stream operator. 
	 * @param os Output stream.
	 * @param _m Matrix. 
	 * @return Output stream reference. 
	 */
	friend std::ostream & operator << (std::ostream & os, const DenseMatrix & _m)
	{
		for (unsigned i = 0; i < _m.m_row; ++i)
		{
			for (unsigned j = 0; j < _m.m_col; ++j)
			{
				os << _m(i,j)<<' ';
			}
			os<<std::endl;
		}
		os<<std::endl;
		return os ;
	}


	/**
	 * Output stream operator. 
	 * @param is Input stream.
	 * @param _m Matrix. 
	 * @return Output stream reference. 
	 */
	friend std::ostream & operator >> (std::ostream & is, const DenseMatrix & _m)
	{
		for (unsigned i = 0; i < _m.m_row; ++i)
		{
			is >> _m[i]; 
		}
		return is;
	}

	/**
	 * Overload operator = (ValueT) for matrix initialization. 
	 * @param _v Value. 
	 * @return Matrix reference. 
	 */
	DenseMatrix & operator = (const ValueT & _v) 
	{
		this->m_idxcount = 0;
		m_pElements[0][0] = _v; 
		return (*this); 
	}

	/**
	 * Overload operator comma  for matrix initialization.
	 * @param _v Value.
	 * @return Matrix reference. 
	 */
	DenseMatrix & operator , (const ValueT & _v)
	{
		++ this->m_idxcount ;
		if (this->m_idxcount >= this->m_row * this->m_col) 
		{
			throw ("Initialization exceeds the matrix dimension bounds."); 
			return (*this);
		}
		unsigned row = m_idxcount / m_col ;
		unsigned col = m_idxcount % m_col ;
		m_pElements[row][col] = _v; 

		return (*this);
	}

	/**
	 * Determinant of 3x3 matrix. Dimension of other matrix is not applicable. The lu_decomposor should be used instead. 
	 */
	ValueT det2 () const
	{
		const 	Vector<Vector<ValueT> > & m = m_pElements ;
		ValueT det = m[0][0] * m[1][1] - m[0][1] * m[1][0]; 
		return det ;
	}

	/**
	 * Determinant of 3x3 matrix. Dimension of other matrix is not applicable. The lu_decomposor should be used instead. 
	 */
	ValueT det3 () const
	{
		const 	Vector<Vector<ValueT> > & m = m_pElements ;
		ValueT det = m[0][0] *(m[2][2]*m[1][1] - m[2][1]* m[1][2]) 
			- m[1][0] * (m[2][2]*m[0][1] - m[2][1]* m[0][2])
			+ m[2][0] * (m[1][2]*m[0][1] - m[1][1]* m[0][2]);
		return det ;
	}

	/**
	 * Inverse matrix of 3x3. 
	 */
	friend DenseMatrix inv3 (const DenseMatrix & m)
	{
		DenseMatrix invm(3,3) ;
        double rdet = 1.0f / m.det3();
        invm[0][0] =   rdet * (m[1][1] * m[2][2] - m[1][2] * m[2][1]) ;
        invm[0][1] = - rdet * (m[0][1] * m[2][2] - m[0][2] * m[2][1]) ;
        invm[0][2] =   rdet * (m[0][1] * m[1][2] - m[0][2] * m[1][1]) ;

        invm[1][0] = - rdet * (m[1][0] * m[2][2] - m[1][2] * m[2][0]) ;
        invm[1][1] =   rdet * (m[0][0] * m[2][2] - m[0][2] * m[2][0]) ;
        invm[1][2] = - rdet * (m[0][0] * m[1][2] - m[0][2] * m[1][0]) ;

        invm[2][0] =   rdet * (m[1][0] * m[2][1] - m[1][1] * m[2][0]) ;
        invm[2][1] = - rdet * (m[0][0] * m[2][1] - m[0][1] * m[2][0]) ;
        invm[2][2] =   rdet * (m[0][0] * m[1][1] - m[0][1] * m[1][0]) ;
		return invm; 
	} 
private:
	void allocate_memory ()
	{
		allocate_memory (m_pElements, m_row, m_col, ValueT(0));
	}
	void allocate_memory (Vector<Vector<ValueT> > & _m, unsigned _r, unsigned _c , const ValueT & _v)
	{
		if (!(_m.size() == 0 ))
		{
			_m.clear();
		}
		
		Vector<ValueT> temp(_c);
		for (unsigned i = 0; i < _r; ++i)
		{
			_m.append(temp);
		}
		for (unsigned i = 0; i < _r ; ++i)
		{
			for (unsigned j = 0; j < _c; ++j)
			{
				_m[i][j] = _v;
			}
		}
	}
	void deallocate_memory()
	{
		//this->deallocate_memory (this->m_pElements, m_row, m_col);
	}

private:
	unsigned m_row, m_col ;

	Vector<Vector<ValueT> > m_pElements ;


	/**
	 * For overloading comma. 
	 */
	unsigned m_idxcount ;  
};

/*
template <typename Value>
struct my_point_map
{
    typedef Vector<Value> Point;
};
*/

#endif
